{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Greatest Common Divisor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import math\n",
    "import random\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gcd_euclid(a, b, debug=False):\n",
    "    i = 1\n",
    "    r0, r1 = a, b\n",
    "    s0, s1 = 1, 0\n",
    "    t0, t1 = 0, 1\n",
    "    print(' ', r0, s0, t0) if debug else None\n",
    "    assert a * s0 + b * t0 == r0\n",
    "    assert a * s1 + b * t1 == r1\n",
    "    while r1 != 0:\n",
    "        q = r0 // r1\n",
    "        print(q, r1, s1, t1) if debug else None\n",
    "        r0, r1 = r1, r0 - q * r1\n",
    "        s0, s1 = s1, s0 - q * s1\n",
    "        t0, t1 = t1, t0 - q * t1\n",
    "        assert r0 > r1\n",
    "        assert a * s0 + b * t0 == r0\n",
    "        assert a * s1 + b * t1 == r1\n",
    "        assert (s1 < 0) == (i % 2 == 0)\n",
    "        assert (t1 < 0) == (i % 2 == 1)\n",
    "        i += 1\n",
    "    print(' ', r1, s1, t1) if debug else None\n",
    "    i -= 1\n",
    "    return (i, r0, s0, t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  126 1 0\n",
      "3 35 0 1\n",
      "1 21 1 -3\n",
      "1 14 -1 4\n",
      "2 7 2 -7\n",
      "  0 -5 18\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(4, 7, 2, -7)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd_euclid(126, 35, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  127 1 0\n",
      "2 45 0 1\n",
      "1 37 1 -2\n",
      "4 8 -1 3\n",
      "1 5 5 -14\n",
      "1 3 -6 17\n",
      "1 2 11 -31\n",
      "2 1 -17 48\n",
      "  0 45 -127\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(7, 1, -17, 48)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd_euclid(127, 45, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "bits = 64\n",
    "A = 26498041357\n",
    "B = 8378459450"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  26498041357 1 0\n",
      "3 8378459450 0 1\n",
      "6 1362663007 1 -3\n",
      "6 202481408 -6 19\n",
      "1 147774559 37 -117\n",
      "2 54706849 -43 136\n",
      "1 38360861 123 -389\n",
      "2 16345988 -166 525\n",
      "2 5668885 455 -1439\n",
      "1 5008218 -1076 3403\n",
      "7 660667 1531 -4842\n",
      "1 383549 -11793 37297\n",
      "1 277118 13324 -42139\n",
      "2 106431 -25117 79436\n",
      "1 64256 63558 -201011\n",
      "1 42175 -88675 280447\n",
      "1 22081 152233 -481458\n",
      "1 20094 -240908 761905\n",
      "10 1987 393141 -1243363\n",
      "8 224 -4172318 13195535\n",
      "1 195 33771685 -106807643\n",
      "6 29 -37944003 120003178\n",
      "1 21 261435703 -826826711\n",
      "2 8 -299379706 946829889\n",
      "1 5 860195115 -2720486489\n",
      "1 3 -1159574821 3667316378\n",
      "1 2 2019769936 -6387802867\n",
      "2 1 -3179344757 10055119245\n",
      "  0 8378459450 -26498041357\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(27, 1, -3179344757, 10055119245)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd_euclid(A, B, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def euclid_qs(a, b):\n",
    "    qs = []\n",
    "    while b != 0:\n",
    "        q = a // b\n",
    "        qs.append(q)\n",
    "        a, b = b, a - q * b\n",
    "    return qs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lehmer_matrix(qs):\n",
    "    q00, q01 = 1, 0\n",
    "    q10, q11 = 0, 1\n",
    "    even = True\n",
    "    for q in qs:\n",
    "        q00, q01, q10, q11 = q10, q11, q00 - q * q10, q01 - q * q11\n",
    "        even = not even\n",
    "    return (q00, q01, q10, q11, even)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = euclid_qs(A, B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3,\n",
       " 6,\n",
       " 6,\n",
       " 1,\n",
       " 2,\n",
       " 1,\n",
       " 2,\n",
       " 2,\n",
       " 1,\n",
       " 7,\n",
       " 1,\n",
       " 1,\n",
       " 2,\n",
       " 1,\n",
       " 1,\n",
       " 1,\n",
       " 1,\n",
       " 10,\n",
       " 8,\n",
       " 1,\n",
       " 6,\n",
       " 1,\n",
       " 2,\n",
       " 1,\n",
       " 1,\n",
       " 1,\n",
       " 2]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-43, 136, 123, -389, False)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lehmer_matrix(Q[0:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(115, -164, -291, 415, True)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lehmer_matrix(Q[5:13])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gcd_small(r0, r1):\n",
    "    assert r0 < 2**bits\n",
    "    assert r1 < 2**bits\n",
    "    q00, q01 = 1, 0\n",
    "    q10, q11 = 0, 1\n",
    "    even = True\n",
    "    while r1 != 0:\n",
    "        q = r0 // r1\n",
    "        r0, r1 = r1, r0 - q * r1\n",
    "        q00, q01, q10, q11 = q10, q11, q * q10 + q00, q * q11 + q01\n",
    "        even = not even\n",
    "    return (q00, q01, q10, q11, even)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lehmer_loop(a, b, q00, q01, q10, q11, even, debug=False):\n",
    "    q0, q1 = None, None\n",
    "    t0, t1 = 0, None\n",
    "    if b != 0:\n",
    "        q0 = a // b\n",
    "        t0 = a - q0 * b\n",
    "    if t0 >= 2**(bits // 2):\n",
    "        while True:\n",
    "            q1 = b // t0\n",
    "            t1 = b - q1 * t0\n",
    "            if t1 < 2**(bits // 2):\n",
    "                # Stopping condition from Cohen et al. algo 10.46 \n",
    "                break\n",
    "            a, b = b, t0\n",
    "            q00, q01, q10, q11 = q10, q11, q0 * q10 + q00, q0 * q11 + q01\n",
    "            print(q0, (q00, q01, q10, q11)) if debug else None\n",
    "            assert q00 < 2**bits and q01 < 2**bits and q11 < 2**bits and q11 < 2**bits\n",
    "            t0, t1 = t1, None\n",
    "            q0, q1 = q1, None\n",
    "            even = not even\n",
    "    return (q00, q01, q10, q11, even)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lehmer_step(A, B, debug=False):\n",
    "    if A < 2**bits:\n",
    "        return gcd_small(A, B)\n",
    "    \n",
    "    s = max(int(math.ceil(math.log2(A))) - bits, 0)\n",
    "    a = A >> s\n",
    "    b = B >> s\n",
    "    print(s, a, b) if debug else None\n",
    "    assert a < 2**bits\n",
    "    assert b < 2**bits\n",
    "    \n",
    "    q00, q01, q10, q11, even = lehmer_loop(a, b, 1, 0, 0, 1, True)\n",
    "    print(q00, q01, q10, q11) if debug else None\n",
    "\n",
    "    # return (q00, q01, q10, q11, even)\n",
    "\n",
    "    if q10 == 0:\n",
    "        # Do a full step: This is done in full precission!\n",
    "        q0 = a // b\n",
    "        q00, q01 = 0, 1\n",
    "        q10, q11 = 1, q0\n",
    "        return (q00, q01, q10, q11, False)\n",
    "\n",
    "    \n",
    "    # Recompute a and b (uses a 2x1 multiply taking the high bits)\n",
    "    s = max(int(math.ceil(math.log2(A))) - 2 * bits, 0)\n",
    "    a = A >> s\n",
    "    b = B >> s\n",
    "    assert a < 2**(2 * bits)\n",
    "    assert b < 2**(2 * bits)\n",
    "    if even:\n",
    "        a, b = q00 * a - q01 * b, q11 * b - q10 * a\n",
    "    else:\n",
    "        a, b = q01 * b - q00 * a, q10 * a - q11 * b\n",
    "    s = max(int(math.ceil(math.log2(a))) - bits, 0)\n",
    "    a >>= s\n",
    "    b >>= s\n",
    "    assert a < 2**bits\n",
    "    assert b < 2**bits\n",
    "    assert a > b\n",
    "\n",
    "    # Itterate once more using new a and b\n",
    "    q00, q01, q10, q11, even = lehmer_loop(a, b, q00, q01, q10, q11, even)\n",
    "\n",
    "    return (q00, q01, q10, q11, even)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6 (1, 0, 0, 1, True)\n",
      "6 (0, 1, 1, -6, False)\n",
      "1 (1, -6, -1, 7, True)\n",
      "2 (-1, 7, 3, -20, False)\n"
     ]
    }
   ],
   "source": [
    "for i in range(2, 6):\n",
    "    print(Q[i - 1], lehmer_matrix(Q[2:i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gcd_lehmer(a, b, debug=False):\n",
    "    r0, r1 = a, b\n",
    "    s0, s1 = 1, 0\n",
    "    t0, t1 = 0, 1\n",
    "    print(r0, s0, t0) if debug else None\n",
    "    assert a * s0 + b * t0 == r0\n",
    "    assert a * s1 + b * t1 == r1\n",
    "    print('t =', t0)\n",
    "    print('t =', t1)\n",
    "    while r1 != 0:\n",
    "        q00, q01, q10, q11, even = lehmer_step(r0, r1)\n",
    "        print(q00, q01, q10, q11, even) if debug else None\n",
    "        if even:\n",
    "            r0, r1 = q00 * r0 - q01 * r1, q11 * r1 - q10 * r0\n",
    "            s0, s1 = q00 * s0 - q01 * s1, q11 * s1 - q10 * s0\n",
    "            t0, t1 = q00 * t0 - q01 * t1, q11 * t1 - q10 * t0\n",
    "        else:\n",
    "            # \n",
    "            r0, r1 = q01 * r1 - q00 * r0, q10 * r0 - q11 * r1\n",
    "            s0, s1 = q01 * s1 - q00 * s0, q10 * s0 - q11 * s1\n",
    "            t0, t1 = q01 * t1 - q00 * t0, q10 * t0 - q11 * t1\n",
    "        assert r0 > r1\n",
    "        assert a * s0 + b * t0 == r0\n",
    "        assert a * s1 + b * t1 == r1\n",
    "        # assert (s1 < 0) == (i % 2 == 0)\n",
    "        # assert (t1 < 0) == (i % 2 == 1)\n",
    "        print('t =', t1)\n",
    "    print(r1, s1, t1) if debug else None\n",
    "    return (i, r0, s0, t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26498041357 1 0\n",
      "t = 0\n",
      "t = 1\n",
      "3179344757 10055119245 8378459450 26498041357 False\n",
      "t = -26498041357\n",
      "0 8378459450 -26498041357\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(5, 1, -3179344757, 10055119245)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd_lehmer(A, B, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0277c215d211565800a282763e0daa0c681edb446debdc8040b335e220ec1f0a\n",
      "0000000000210149c777d4845be3b922780375a9f0abda3e085b78a3e366a8fa\n",
      "0000000000000000000000000029be0ec7dc762158e31f89fd4f5e834390ec8e\n",
      "000000000000000000000000000000000000000000000000455b2a337ec2b649\n",
      "00000000000000000000000000000000000000052f905beba02d4e4bf1eb761e\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "g = random.randint(0, 2**150)\n",
    "a = random.randint(0, 2**100) * g\n",
    "b = random.randint(0, 2**64) * g\n",
    "a, b = max(a, b), min(a, b)\n",
    "print(a.to_bytes(32, 'big').hex())\n",
    "print(b.to_bytes(32, 'big').hex())\n",
    "i, g, s, t  = gcd_euclid(a, b)\n",
    "print(g.to_bytes(32, 'big').hex())\n",
    "print(abs(s).to_bytes(32, 'big').hex())\n",
    "print(abs(t).to_bytes(32, 'big').hex())\n",
    "print(i % 2 == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way of looking at Lehmer's GCD is as a process of accumulating Euclid quotients into matrices. We could greedily accumulate coefficients untill we overflow the word size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lehmer_accumulate(qs):\n",
    "    q00, q01, q10, q11 = 1, 0, 0, 1\n",
    "    qss = []\n",
    "    for q in qs:\n",
    "        q20 = q00 + q * q10\n",
    "        q21 = q01 + q * q11\n",
    "        if q20 >= 2**64 or q21 >= 2**64:\n",
    "            print(qss)\n",
    "            yield((q00, q01, q10, q11))\n",
    "            q00, q01, q10, q11 = 1, 0, 0, 1\n",
    "            q20 = q00 + q * q10\n",
    "            q21 = q01 + q * q11\n",
    "            qss = [q]\n",
    "        q00, q10 = q10, q20\n",
    "        q01, q11 = q11, q21\n",
    "        qss += [q]\n",
    "    print(qss)\n",
    "    yield (q00, q01, q10, q11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this assumes the quotients themselves always fit a word, which is not necessarily true."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 344,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = 0xc3e61c8ad68eb12ba57d7e9723e3def765de4102d10e6ce8372bae1a88a51115\n",
      "B = 0x6935d680eb6f04f9c9b09a1501baee776267db67245c4aff2d551542eba3247b\n"
     ]
    }
   ],
   "source": [
    "A = random.randint(0, 2**256)\n",
    "B = random.randint(0, A)\n",
    "print('A = 0x' + A.to_bytes(32, 'big').hex())\n",
    "print('B = 0x' + B.to_bytes(32, 'big').hex())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 345,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = 0xfea5a792d0a17b24827908e5524bcceec3ec6a92a7a42eac3b93e2bb351cf4f2\n",
    "B = 0x00028735553c6c798ed1ffb8b694f8f37b672b1bab7f80c4e6f4c0e710c79fb4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 346,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'25785 2 3 1 1 1 29 23 4 2 2 2 3 1 1 1 6 2 5 1 3 1 3 2 4 2 3 5 1 4 1 1 4 2 4 7 3 3 1 3 3 2 2 11 1 39 1 1 2 1 2 2 1 22 1 7 1 1 1 1 1 1 3 6 1 3 1 1 310 1 1 2 1 1 1 2 2 1 1 8 2 10 1 30 21 1 1 6 4 1 1 3 1 2 1 5 1 2 1 1 1 1 1 31 1 2 1 1 1 6 1 7 4 1 14 1 6 3 1 6 1 2 1 1 1 4 1 6 1 2 2 1 2 1 26 1 3 4 2 1 2 2 2 3 4 1 1 1 39 16 1 1 2 7 1 64 1 1 28'"
      ]
     },
     "execution_count": 346,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "' '.join([str(x) for x in list(euclid_qs(A, B))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 340,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "159\n",
      "g = 0x0000000000000000000000000000000000000000000000000000000000000002\n",
      "u = 0x00000b5a5ecb4dfc4ea08773d0593986592959a646b2f97655ed839928274ebb\n",
      "v = 0x0477865490d3994853934bf7eae7dad9afac55ccbf412a60c18fc9bea58ec8ba\n"
     ]
    }
   ],
   "source": [
    "(c, g, u, v) = gcd_euclid(A, B)\n",
    "print(c)\n",
    "print('g = 0x' + g.to_bytes(32, 'big').hex())\n",
    "print('u = 0x' + abs(u).to_bytes(32, 'big').hex())\n",
    "print('v = 0x' + abs(v).to_bytes(32, 'big').hex())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 343,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'25785 2 3 1 1 1 29 23 4 1 2 2 2 3 1 1 1 6 2 5 1 3 1 3 2 4 2 3 5 1 4 1 1'"
      ]
     },
     "execution_count": 343,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qs = [25785, 2, 3, 1, 1, 1, 29, 23, 4, 1]\n",
    "qs += [2, 2, 2, 3, 1,\n",
    "1,\n",
    "1,\n",
    "6,\n",
    "2,\n",
    "5,\n",
    "1,\n",
    "3,\n",
    "1,\n",
    "3,\n",
    "2,\n",
    "4,\n",
    "2,\n",
    "3,\n",
    "5,\n",
    "1,\n",
    "4,\n",
    "1,\n",
    "1]\n",
    "' '.join([str(x) for x in list(qs)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "euclid_qs(A, B) == qs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 347,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[25785, 2, 3, 1, 1, 1, 29, 23, 4]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[(17068, 440105889, 69013, 1779530567)]"
      ]
     },
     "execution_count": 347,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(lehmer_accumulate([25785, 2, 3, 1, 1, 1, 29, 23, 4]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 348,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2, 2, 2, 3, 1, 1, 1, 6, 2, 5, 1, 3, 1, 3, 2, 4, 2, 3, 5, 1, 4, 1, 1]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[(293461651, 707205579, 536681752, 1293335357)]"
      ]
     },
     "execution_count": 348,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(lehmer_accumulate([2, 2, 2, 3, 1,\n",
    "1,\n",
    "1,\n",
    "6,\n",
    "2,\n",
    "5,\n",
    "1,\n",
    "3,\n",
    "1,\n",
    "3,\n",
    "2,\n",
    "4,\n",
    "2,\n",
    "3,\n",
    "5,\n",
    "1,\n",
    "4,\n",
    "1,\n",
    "1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 349,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[25785, 2, 3, 1, 1, 1, 29, 23, 4, 2, 2, 2, 3, 1, 1, 1, 6, 2, 5, 1, 3, 1, 3, 2, 4, 2, 3, 5, 1, 4, 1, 1]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[(53815182082795, 1387648145784196032, 98417037135777, 2537726600737394947)]"
      ]
     },
     "execution_count": 349,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(lehmer_accumulate([25785, 2, 3, 1, 1, 1, 29, 23, 4, 2, 2, 2, 3, 1,\n",
    "1,\n",
    "1,\n",
    "6,\n",
    "2,\n",
    "5,\n",
    "1,\n",
    "3,\n",
    "1,\n",
    "3,\n",
    "2,\n",
    "4,\n",
    "2,\n",
    "3,\n",
    "5,\n",
    "1,\n",
    "4,\n",
    "1,\n",
    "1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 353,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "53815182082795"
      ]
     },
     "execution_count": 353,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "17068 * 293461651 + 69013 * 707205579"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[25785, 2, 3, 1, 1, 1, 29, 23, 4, 2, 2, 2, 3, 1, 1, 1, 6, 2, 5, 1, 3, 1, 3, 2, 4, 2, 3, 5, 1, 4, 1, 1, 4]\n",
      "(98417037135777, 2537726600737394947, 447483330625903, 11538554548733775820)\n",
      "[2, 2, 4, 7, 3, 3, 1, 3, 3, 2, 2, 11, 1, 39, 1, 1, 2, 1, 2, 2, 1, 22, 1, 7, 1, 1, 1, 1, 1, 1, 3, 6, 1, 3, 1, 1, 310, 1, 1, 2, 1, 1, 1]\n",
      "(1977056207563995501, 4432023614448106363, 3130294868778833684, 7017271803166199793)\n",
      "[2, 2, 2, 1, 1, 8, 2, 10, 1, 30, 21, 1, 1, 6, 4, 1, 1, 3, 1, 2, 1, 5, 1, 2, 1, 1, 1, 1, 1, 31, 1, 2, 1, 1, 1, 6, 1, 7, 4, 1, 14, 1, 6]\n",
      "(470618561295672659, 1127410926475022632, 3264553286537374601, 7820543743895941587)\n",
      "[3, 3, 1, 6, 1, 2, 1, 1, 1, 4, 1, 6, 1, 2, 2, 1, 2, 1, 26, 1, 3, 4, 2, 1, 2, 2, 2, 3, 4, 1, 1, 1, 39, 16, 1, 1, 2, 7, 1, 64, 1, 1]\n",
      "(288562567703758841, 1116899625196145864, 572745072480815301, 2216845940473758943)\n",
      "[28, 28]\n",
      "(0, 1, 1, 28)\n"
     ]
    }
   ],
   "source": [
    "for q in lehmer_accumulate(euclid_qs(A, B)):\n",
    "    print(q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lehmer_loop(r0, r1, debug=False):\n",
    "    assert r0 < 2**bits\n",
    "    assert r1 < 2**bits\n",
    "    assert r0 >= 2**(bits - 1)\n",
    "    assert r0 >= r1\n",
    "    a = [r0, r1]\n",
    "    q = [None]\n",
    "    u = [1, 0]\n",
    "    v = [0, 1]\n",
    "    while True:\n",
    "        q += [a[-2] // a[-1]]\n",
    "        a += [a[-2] - q[-1] * a[-1]]\n",
    "        u += [u[-2] + q[-1] * u[-1]]\n",
    "        v += [v[-2] + q[-1] * v[-1]]\n",
    "        # Stopping condition\n",
    "        if a[-1] < 2**(bits // 2):\n",
    "            break\n",
    "            \n",
    "    print(a)\n",
    "    print(q)\n",
    "    print(u)\n",
    "    print(v)\n",
    "        \n",
    "    i = len(a) - 4\n",
    "    assert a[i + 2] >= 2**(bits // 2)\n",
    "    assert a[i + 3] < 2**(bits // 2)\n",
    "    print('i =', i)\n",
    "    k = i\n",
    "    \n",
    "    print([a[i], a[i+1], a[i+2], a[i+3]])\n",
    "    print([u[i], u[i+1], u[i+2], u[i+3]])\n",
    "    print([v[i], v[i+1], v[i+2], v[i+3]])\n",
    "    \n",
    "    # Test right hand side for i + 1\n",
    "    if i % 2 == 0:\n",
    "        # Test i + 1 (odd)\n",
    "        assert a[i + 2] >= v[i + 2]\n",
    "        if a[i + 1] - a[i + 2] >= u[i + 2] - u[i + 1]:\n",
    "            # Test i + 3 (odd)\n",
    "            if a[i + 3] >= u[i + 3] and a[i + 2] - a[i + 3] >= u[i + 3] - u[i + 2]:\n",
    "                print('k = i + 2')\n",
    "                k = i + 2\n",
    "            else:\n",
    "                print('k = i + 1')\n",
    "                k = i + 1\n",
    "        else:\n",
    "            k = i\n",
    "            print('k = i')\n",
    "    else:\n",
    "        # Test i + 1 (even)\n",
    "        assert a[i + 2] >= u[i + 2]\n",
    "        if a[i + 1] - a[i + 2] >= v[i + 2] - v[i + 1]:\n",
    "            # Test i + 2 (odd)\n",
    "            if a[i + 3] >= v[i + 3] and a[i + 2] - a[i + 3] >= v[i + 3] - v[i + 2]:\n",
    "                k = i + 2\n",
    "                print('k = i + 2')\n",
    "            else:\n",
    "                k = i + 1\n",
    "                print('k = i + 1')\n",
    "        else:\n",
    "            k = i\n",
    "            print('k = i')\n",
    "    print(k)\n",
    "    print(q[k])\n",
    "    \n",
    "    print(a[:k+2])\n",
    "    print(q[:k+1])\n",
    "    print(u[:k+2])\n",
    "    print(v[:k+2])\n",
    "    return (u[k], u[k+1], v[k], v[k+1], k % 2 == 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 273,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[12664002055393416730, 3174291274789560733, 3141128231024734531, 33163043764826202, 23802117131071543, 9360926633754659, 5080263863562225, 4280662770192434, 799601093369791, 282657303343479, 234286486682833, 48370816660646, 40803220040249, 7567596620397, 2965236938264, 1637122743869, 1328114194395, 309008549474, 92079996499, 32768559977, 26542876545, 6225683432, 1640142817]\n",
      "[None, 3, 1, 94, 1, 2, 1, 1, 5, 2, 1, 4, 1, 5, 2, 1, 1, 4, 3, 2, 1, 4]\n",
      "[1, 0, 1, 1, 95, 96, 287, 383, 670, 3733, 8136, 11869, 55612, 67481, 393017, 853515, 1246532, 2100047, 9646720, 31040207, 71727134, 102767341, 482796498]\n",
      "[0, 1, 3, 4, 379, 383, 1145, 1528, 2673, 14893, 32459, 47352, 221867, 269219, 1567962, 3405143, 4973105, 8378248, 38486097, 123836539, 286159175, 409995714, 1926142031]\n",
      "i = 19\n",
      "[32768559977, 26542876545, 6225683432, 1640142817]\n",
      "[31040207, 71727134, 102767341, 482796498]\n",
      "[123836539, 286159175, 409995714, 1926142031]\n",
      "k = i + 1\n",
      "20\n",
      "1\n",
      "[12664002055393416730, 3174291274789560733, 3141128231024734531, 33163043764826202, 23802117131071543, 9360926633754659, 5080263863562225, 4280662770192434, 799601093369791, 282657303343479, 234286486682833, 48370816660646, 40803220040249, 7567596620397, 2965236938264, 1637122743869, 1328114194395, 309008549474, 92079996499, 32768559977, 26542876545, 6225683432]\n",
      "[None, 3, 1, 94, 1, 2, 1, 1, 5, 2, 1, 4, 1, 5, 2, 1, 1, 4, 3, 2, 1]\n",
      "[1, 0, 1, 1, 95, 96, 287, 383, 670, 3733, 8136, 11869, 55612, 67481, 393017, 853515, 1246532, 2100047, 9646720, 31040207, 71727134, 102767341]\n",
      "[0, 1, 3, 4, 379, 383, 1145, 1528, 2673, 14893, 32459, 47352, 221867, 269219, 1567962, 3405143, 4973105, 8378248, 38486097, 123836539, 286159175, 409995714]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(71727134, 102767341, 286159175, 409995714, True)"
      ]
     },
     "execution_count": 273,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bits = 64\n",
    "s = max(int(math.ceil(math.log2(A))) - bits, 0)\n",
    "a = A >> s\n",
    "b = B >> s\n",
    "lehmer_loop(a, b, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 334,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lehmer_loop_3(a0, a1, debug=False):\n",
    "    assert a0 < 2**bits\n",
    "    assert a1 < 2**bits\n",
    "    assert a0 >= 2**(bits - 1)\n",
    "    assert a0 >= a1\n",
    "    LIMIT = 2**(bits // 2)\n",
    "    u0, u1 = 1, 0\n",
    "    v0, v1 = 0, 1\n",
    "    even = True\n",
    "    \n",
    "    # Compute a2\n",
    "    q = a0 // a1\n",
    "    a2 = a0 - q * a1\n",
    "    u2 = u0 + q * u1\n",
    "    v2 = v0 + q * v1\n",
    "    if a2 < LIMIT:\n",
    "        return (u0, u1, v0, v1, True)\n",
    "    \n",
    "    # Compute a3\n",
    "    q = a1 // a2\n",
    "    a3 = a1 - q * a2\n",
    "    u3 = u1 + q * u2\n",
    "    v3 = v1 + q * v2\n",
    "    \n",
    "    # Itterate\n",
    "    while a3 >= LIMIT:\n",
    "        q = a2 // a3\n",
    "        a0, a1, a2, a3 = a1, a2, a3, a2 - q * a3\n",
    "        u0, u1, u2, u3 = u1, u2, u3, u2 + q * u3\n",
    "        v0, v1, v2, v3 = v1, v2, v3, v2 + q * v3\n",
    "        even = !even\n",
    "    \n",
    "    assert a2 >= LIMIT\n",
    "    assert a3 < LIMIT\n",
    "    \n",
    "    print([a0, a1, a2, a3]) if debug else None\n",
    "    print([u0, u1, u2, u3]) if debug else None\n",
    "    print([v0, v1, v2, v3]) if debug else None\n",
    "    \n",
    "    i = 0\n",
    "    if even:\n",
    "        # Test i + 1 (odd)\n",
    "        assert a2 >= v2\n",
    "        if a1 - a2 >= u2 + u1:\n",
    "            # Test i + 3 (odd)\n",
    "            if a3 >= u3 and a2 - a3 >= u3 + u2:\n",
    "                print('k = i + 2') if debug else None\n",
    "                k = i + 2\n",
    "                return 2\n",
    "                return (u2, u3, v2, v3, True)\n",
    "            else:\n",
    "                print('k = i + 1') if debug else None\n",
    "                k = i + 1\n",
    "                return 1\n",
    "                return (u1, u2, v1, v2, False)\n",
    "        else:\n",
    "            k = i\n",
    "            print('k = i') if debug else None\n",
    "            return 0\n",
    "            return (u0, u1, v0, v1, True)\n",
    "    else:\n",
    "        # Test i + 1 (even)\n",
    "        assert a2 >= u2\n",
    "        if a1 - a2 >= v2 + v1:\n",
    "            # Test i + 2 (odd)\n",
    "            if a3 >= v3 and a2 - a3 >= v3 + v2:\n",
    "                k = i + 2\n",
    "                print('k = i + 2') if debug else None\n",
    "                return 2\n",
    "                return (u2, u3, v2, v3, False)\n",
    "            else:\n",
    "                k = i + 1\n",
    "                print('k = i + 1') if debug else None\n",
    "                return 1\n",
    "                return (u1, u2, v1, v2, True)\n",
    "        else:\n",
    "            k = i\n",
    "            print('k = i') if debug else None\n",
    "            return 0\n",
    "            return (u0, u1, v0, v1, False)\n",
    "    print(k)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 335,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[59590115859, 7697006119, 5711073026, 1985933093]\n",
      "[29647479, 67161500, 499777979, 566939479]\n",
      "[87622597, 198494619, 1477084930, 1675579549]\n",
      "k = i + 2\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 335,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bits = 64\n",
    "s = max(int(math.ceil(math.log2(A))) - bits, 0)\n",
    "a = A >> s\n",
    "b = B >> s\n",
    "lehmer_loop_3(a, b, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 336,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 85., 152., 763.]), array([0, 1, 2, 3]), <a list of 3 Patch objects>)"
      ]
     },
     "execution_count": 336,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEi5JREFUeJzt3V2MXdd53vH/E1GyHTsV9TFhCJIOFZhIoAaxzAwEGi4C16wLSQlMAVUEGUVECywmaNXGjgs0bC4apOiFDRRRo7ZQQJhuqcBxpCh2xTpKGoJSEPRCTEayLOvDrsZqFJKgxLEs0XGUL6ZvL86iPWKGPPtwZjic5f8PODhrr7322e/SoZ7Zs+Z8pKqQJPXre1a7AEnSyjLoJalzBr0kdc6gl6TOGfSS1DmDXpI6Nyjok/x8kmeTPJPks0nemuS6JEeSzCV5IMkVbexb2vZc2791JScgSTq/sUGfZBPwc8B0Vf0ocBlwB/BJ4J6qehfwGrCnHbIHeK3139PGSZJWydClm3XA25KsA74XOAF8AHio7T8A3Nrau9o2bf/OJFmeciVJk1o3bkBVHU/yH4E/Bf4C+H3gCeD1qjrdhh0DNrX2JuBoO/Z0klPANcDXz3WOa6+9trZu3Xqhc5Ck70pPPPHE16tqaty4sUGf5CpGV+nXAa8DvwXctNQCk8wAMwDvfOc7mZ2dXepDStJ3lSQvDRk3ZOnmHwH/t6rmq+pvgM8B7wPWt6UcgM3A8dY+DmxpRawDrgRePftBq2pfVU1X1fTU1NgfSJKkCzQk6P8U2JHke9ta+07gOeAx4LY2ZjfwcGsfbNu0/Y+Wn5wmSatmbNBX1RFGf1R9EvhyO2Yf8AvAx5PMMVqD398O2Q9c0/o/DuxdgbolSQPlUrjYnp6eLtfoJWkySZ6oqulx43xnrCR1zqCXpM4Z9JLUOYNekjpn0EtS58a+M1ZSH7bu/Z3VLkGL+JNP/OSKn8MreknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1bmzQJ/nhJE8tuH0zyceSXJ3kUJIX2v1VbXyS3JtkLsnTSbav/DQkSecy5MvBv1pVN1TVDcCPA28An2f0pd+Hq2obcJjvfAn4zcC2dpsB7luJwiVJw0y6dLMT+FpVvQTsAg60/gPAra29C7i/Rh4H1ifZuCzVSpImNmnQ3wF8trU3VNWJ1n4Z2NDam4CjC4451vokSatgcNAnuQL4EPBbZ++rqgJqkhMnmUkym2R2fn5+kkMlSROY5Ir+ZuDJqnqlbb9yZkmm3Z9s/ceBLQuO29z63qSq9lXVdFVNT01NTV65JGmQSYL+w3xn2QbgILC7tXcDDy/ov7O9+mYHcGrBEo8k6SIb9J2xSd4OfBD42QXdnwAeTLIHeAm4vfU/AtwCzDF6hc5dy1atJGlig4K+qv4cuOasvlcZvQrn7LEF3L0s1UmSlsx3xkpS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6tygoE+yPslDSb6S5Pkk701ydZJDSV5o91e1sUlyb5K5JE8n2b6yU5Aknc/QK/pfBX6vqn4EeDfwPLAXOFxV24DDbRvgZmBbu80A9y1rxZKkiYwN+iRXAj8B7Aeoqr+uqteBXcCBNuwAcGtr7wLur5HHgfVJNi575ZKkQYZc0V8HzAP/LckXk3wqyduBDVV1oo15GdjQ2puAowuOP9b6JEmrYEjQrwO2A/dV1XuAP+c7yzQAVFUBNcmJk8wkmU0yOz8/P8mhkqQJDAn6Y8CxqjrSth9iFPyvnFmSafcn2/7jwJYFx29ufW9SVfuqarqqpqempi60fknSGGODvqpeBo4m+eHWtRN4DjgI7G59u4GHW/sgcGd79c0O4NSCJR5J0kW2buC4fwV8JskVwIvAXYx+SDyYZA/wEnB7G/sIcAswB7zRxkqSVsmgoK+qp4DpRXbtXGRsAXcvsS5J0jLxnbGS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wYFfZI/SfLlJE8lmW19Vyc5lOSFdn9V60+Se5PMJXk6yfaVnIAk6fwmuaL/h1V1Q1Wd+ZLwvcDhqtoGHG7bADcD29ptBrhvuYqVJE1uKUs3u4ADrX0AuHVB//018jiwPsnGJZxHkrQEQ4O+gN9P8kSSmda3oapOtPbLwIbW3gQcXXDssdb3JklmkswmmZ2fn7+A0iVJQ6wbOO4fVNXxJN8PHErylYU7q6qS1CQnrqp9wD6A6enpiY6VJA036Iq+qo63+5PA54EbgVfOLMm0+5Nt+HFgy4LDN7c+SdIqGBv0Sd6e5PvOtIF/DDwDHAR2t2G7gYdb+yBwZ3v1zQ7g1IIlHknSRTZk6WYD8PkkZ8b/RlX9XpI/Bh5Msgd4Cbi9jX8EuAWYA94A7lr2qiVJg40N+qp6EXj3Iv2vAjsX6S/g7mWpTpK0ZL4zVpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wYHfZLLknwxyRfa9nVJjiSZS/JAkita/1va9lzbv3VlSpckDTHJFf1HgecXbH8SuKeq3gW8Buxp/XuA11r/PW2cJGmVDAr6JJuBnwQ+1bYDfAB4qA05ANza2rvaNm3/zjZekrQKhl7R/yfg3wD/r21fA7xeVafb9jFgU2tvAo4CtP2n2vg3STKTZDbJ7Pz8/AWWL0kaZ2zQJ/kp4GRVPbGcJ66qfVU1XVXTU1NTy/nQkqQF1g0Y8z7gQ0luAd4K/D3gV4H1Sda1q/bNwPE2/jiwBTiWZB1wJfDqslcuSRpk7BV9Vf3bqtpcVVuBO4BHq+qfAo8Bt7Vhu4GHW/tg26btf7SqalmrliQNtpTX0f8C8PEkc4zW4Pe3/v3ANa3/48DepZUoSVqKIUs331ZVfwD8QWu/CNy4yJi/BH56GWqTJC0D3xkrSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzY4M+yVuT/FGSLyV5Nskvt/7rkhxJMpfkgSRXtP63tO25tn/ryk5BknQ+Q67o/wr4QFW9G7gBuCnJDuCTwD1V9S7gNWBPG78HeK3139PGSZJWydigr5Fvtc3L262ADwAPtf4DwK2tvatt0/bvTJJlq1iSNJFBa/RJLkvyFHASOAR8DXi9qk63IceATa29CTgK0PafAq5Z5DFnkswmmZ2fn1/aLCRJ5zQo6Kvqb6vqBmAzcCPwI0s9cVXtq6rpqpqemppa6sNJks5holfdVNXrwGPAe4H1Sda1XZuB4619HNgC0PZfCby6LNVKkiY25FU3U0nWt/bbgA8CzzMK/NvasN3Aw619sG3T9j9aVbWcRUuShls3fggbgQNJLmP0g+HBqvpCkueA30zyH4AvAvvb+P3AryeZA74B3LECdUuSBhob9FX1NPCeRfpfZLRef3b/XwI/vSzVSZKWzHfGSlLnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknq3JAvB9+S5LEkzyV5NslHW//VSQ4leaHdX9X6k+TeJHNJnk6yfaUnIUk6tyFX9KeBf11V1wM7gLuTXA/sBQ5X1TbgcNsGuBnY1m4zwH3LXrUkabCxQV9VJ6rqydb+M+B5YBOwCzjQhh0Abm3tXcD9NfI4sD7JxmWvXJI0yERr9Em2Au8BjgAbqupE2/UysKG1NwFHFxx2rPVJklbB4KBP8g7gt4GPVdU3F+6rqgJqkhMnmUkym2R2fn5+kkMlSRMYFPRJLmcU8p+pqs+17lfOLMm0+5Ot/ziwZcHhm1vfm1TVvqqarqrpqampC61fkjTGkFfdBNgPPF9Vv7Jg10Fgd2vvBh5e0H9ne/XNDuDUgiUeSdJFtm7AmPcBPwN8OclTre8XgU8ADybZA7wE3N72PQLcAswBbwB3LWvFkqSJjA36qvrfQM6xe+ci4wu4e4l1SZKWie+MlaTOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHVuyDtjpYlt3fs7q12CpMYreknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6tyQLwf/dJKTSZ5Z0Hd1kkNJXmj3V7X+JLk3yVySp5NsX8niJUnjDbmi/+/ATWf17QUOV9U24HDbBrgZ2NZuM8B9y1OmJOlCjQ36qvpD4Btnde8CDrT2AeDWBf3318jjwPokG5erWEnS5C50jX5DVZ1o7ZeBDa29CTi6YNyx1idJWiVL/mNsVRVQkx6XZCbJbJLZ+fn5pZYhSTqHCw36V84sybT7k63/OLBlwbjNre/vqKp9VTVdVdNTU1MXWIYkaZwLDfqDwO7W3g08vKD/zvbqmx3AqQVLPJKkVTD2i0eSfBZ4P3BtkmPALwGfAB5Msgd4Cbi9DX8EuAWYA94A7lqBmiVJExgb9FX14XPs2rnI2ALuXmpRk/CbjCTp/HxnrCR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzq1I0Ce5KclXk8wl2bsS55AkDbPsQZ/kMuC/AjcD1wMfTnL9cp9HkjTMSlzR3wjMVdWLVfXXwG8Cu1bgPJKkAVYi6DcBRxdsH2t9kqRVsG61TpxkBphpm99K8tULfKhrga8vT1WrzrlcenqZBziXS1I+uaS5/OCQQSsR9MeBLQu2N7e+N6mqfcC+pZ4syWxVTS/1cS4FzuXS08s8wLlcqi7GXFZi6eaPgW1JrktyBXAHcHAFziNJGmDZr+ir6nSSfwn8L+Ay4NNV9exyn0eSNMyKrNFX1SPAIyvx2ItY8vLPJcS5XHp6mQc4l0vVis8lVbXS55AkrSI/AkGSOrdmgn7cxyokeUuSB9r+I0m2Xvwqhxkwl48kmU/yVLv9s9Woc5wkn05yMskz59ifJPe2eT6dZPvFrnGoAXN5f5JTC56Tf3exaxwiyZYkjyV5LsmzST66yJg18bwMnMtaeV7emuSPknypzeWXFxmzchlWVZf8jdEfdb8G/BBwBfAl4PqzxvwL4Nda+w7ggdWuewlz+QjwX1a71gFz+QlgO/DMOfbfAvwuEGAHcGS1a17CXN4PfGG16xwwj43A9tb+PuD/LPLva008LwPnslaelwDvaO3LgSPAjrPGrFiGrZUr+iEfq7ALONDaDwE7k+Qi1jhUNx8RUVV/CHzjPEN2AffXyOPA+iQbL051kxkwlzWhqk5U1ZOt/WfA8/zdd6aviedl4FzWhPbf+ltt8/J2O/sPpCuWYWsl6Id8rMK3x1TVaeAUcM1FqW4yQz8i4p+0X6sfSrJlkf1rQW8fh/He9qv37yb5+6tdzDjtV//3MLp6XGjNPS/nmQuskeclyWVJngJOAoeq6pzPy3Jn2FoJ+u82/xPYWlU/BhziOz/ltXqeBH6wqt4N/Gfgf6xyPeeV5B3AbwMfq6pvrnY9SzFmLmvmeamqv62qGxh9WsCNSX70Yp17rQT9kI9V+PaYJOuAK4FXL0p1kxk7l6p6tar+qm1+Cvjxi1Tbchv0cRhrQVV988yv3jV6n8jlSa5d5bIWleRyRsH4mar63CJD1szzMm4ua+l5OaOqXgceA246a9eKZdhaCfohH6twENjd2rcBj1b7q8YlZuxczlov/RCjtcm16CBwZ3uVxw7gVFWdWO2iLkSSHzizXprkRkb/71xyFxKtxv3A81X1K+cYtiaelyFzWUPPy1SS9a39NuCDwFfOGrZiGbZqn145iTrHxyok+ffAbFUdZPQP4teTzDH6o9odq1fxuQ2cy88l+RBwmtFcPrJqBZ9Hks8yetXDtUmOAb/E6I9MVNWvMXp39C3AHPAGcNfqVDregLncBvzzJKeBvwDuuEQvJN4H/Azw5bYeDPCLwDthzT0vQ+ayVp6XjcCBjL6Y6XuAB6vqCxcrw3xnrCR1bq0s3UiSLpBBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5/4/FCAwiSajGfIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bits = 64\n",
    "ix = []\n",
    "for i in range(1000):\n",
    "    A = random.randint(0, 2**256)\n",
    "    B = random.randint(0, A)\n",
    "\n",
    "    s = max(int(math.ceil(math.log2(A))) - bits, 0)\n",
    "    a = A >> s\n",
    "    b = B >> s\n",
    "    ix += [lehmer_loop_3(a, b)]\n",
    "plt.hist(ix, [0, 1, 2, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hypothesis.** Any 256-bit GCD can be computed using at most 5 Lehmer matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 309,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 4.,  0.,  0.,  0.,  0., 13.,  0.,  0.,  0., 83.]),\n",
       " array([0. , 0.2, 0.4, 0.6, 0.8, 1. , 1.2, 1.4, 1.6, 1.8, 2. ]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 309,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEZJJREFUeJzt3X+MHGd9x/H3pzHhR0DEIYfrJoCNiEChKgFOEb9EIYESSMGuilAQrQx15dIChVK1mCKVtqrURKoKVK1aWYTWSDQkhNCk/CquCUItjeESAvlFiGMSiOXER0iAgASEfvvHjmExd9nZu9098/B+SaedeeaZm2+enXxuPLOzk6pCkvSz7xfWugBJ0mQY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGrJvlxk4++eTatGnTLDcpST/zrr766q9X1dyofjMN9E2bNrGwsDDLTUrSz7wkt/fp5ykXSWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqxEzvFJWktbRp50fWZLu3nX/uTLbjEbokNcJAl6RGGOiS1IhegZ7kj5LckOT6JBcleUiSzUn2Jdmf5OIkx0+7WEnS8kYGepJTgD8E5qvql4HjgPOAC4B3VNUTgHuA7dMsVJL0wPqeclkHPDTJOuBhwCHgLODSbvluYOvky5Mk9TUy0KvqIPC3wFcZBPk3gauBe6vq/q7bHcApS62fZEeShSQLi4uLk6lakvRT+pxyWQ9sATYDvwScAJzTdwNVtauq5qtqfm5u5BOUJEkr1OeUywuAr1TVYlX9ALgMeDZwYncKBuBU4OCUapQk9dAn0L8KPCPJw5IEOBu4EbgSeHnXZxtw+XRKlCT10ecc+j4GFz+vAa7r1tkFvAV4c5L9wKOAC6dYpyRphF7f5VJVbwfeflTzAeDMiVckSVoR7xSVpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEX2eKfrEJNcO/XwryZuSnJRkT5Jbutf1syhYkrS0Pk8surmqzqiqM4CnA98FPgTsBPZW1WnA3m5ekrRGxj3lcjZwa1XdDmwBdnftu4GtkyxMkjSecQP9POCibnpDVR3qpu8ENkysKknS2HoHepLjgZcBHzh6WVUVUMustyPJQpKFxcXFFRcqSXpg4xyhvxi4pqru6ubvSrIRoHs9vNRKVbWrquaran5ubm511UqSljVOoL+SH59uAbgC2NZNbwMun1RRkqTx9Qr0JCcALwQuG2o+H3hhkluAF3TzkqQ1sq5Pp6r6DvCoo9ruZvCpF0nSMcA7RSWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5Jjej7xKITk1ya5EtJbkryzCQnJdmT5Jbudf20i5UkLa/vEfq7gI9X1ZOApwA3ATuBvVV1GrC3m5ckrZGRgZ7kkcBzgQsBqur7VXUvsAXY3XXbDWydVpGSpNH6HKFvBhaBf0ny+STv7h4avaGqDnV97gQ2TKtISdJofQJ9HfA04J+q6qnAdzjq9EpVFVBLrZxkR5KFJAuLi4urrVeStIw+gX4HcEdV7evmL2UQ8Hcl2QjQvR5eauWq2lVV81U1Pzc3N4maJUlLGBnoVXUn8LUkT+yazgZuBK4AtnVt24DLp1KhJKmXdT37vQF4X5LjgQPAaxj8MbgkyXbgduAV0ylRktRHr0CvqmuB+SUWnT3ZciRJK+WdopLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRvR6wEWS24BvAz8E7q+q+SQnARcDm4DbgFdU1T3TKVOSNMo4R+jPr6ozqurIk4t2Anur6jRgbzcvSVojqznlsgXY3U3vBrauvhxJ0kr1DfQCPpHk6iQ7urYNVXWom74T2DDx6iRJvfU6hw48p6oOJnk0sCfJl4YXVlUlqaVW7P4A7AB47GMfu6piJUnL63WEXlUHu9fDwIeAM4G7kmwE6F4PL7Purqqar6r5ubm5yVQtSfopIwM9yQlJHnFkGvg14HrgCmBb120bcPm0ipQkjdbnlMsG4ENJjvT/t6r6eJLPAZck2Q7cDrxiemVKkkYZGehVdQB4yhLtdwNnT6MoSdL4vFNUkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktSI3oGe5Lgkn0/y4W5+c5J9SfYnuTjJ8dMrU5I0yjhH6G8EbhqavwB4R1U9AbgH2D7JwiRJ4+kV6ElOBc4F3t3NBzgLuLTrshvYOo0CJUn99D1Cfyfwp8D/dfOPAu6tqvu7+TuAUyZcmyRpDCMDPcmvA4er6uqVbCDJjiQLSRYWFxdX8iskST30OUJ/NvCyJLcB72dwquVdwIlJ1nV9TgUOLrVyVe2qqvmqmp+bm5tAyZKkpYwM9Kp6a1WdWlWbgPOAT1bVq4ArgZd33bYBl0+tSknSSKv5HPpbgDcn2c/gnPqFkylJkrQS60Z3+bGq+hTwqW76AHDm5EuSJK2Ed4pKUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhrR5yHRD0ny2SRfSHJDkr/s2jcn2Zdkf5KLkxw//XIlScvpc4T+PeCsqnoKcAZwTpJnABcA76iqJwD3ANunV6YkaZQ+D4muqrqvm31Q91PAWcClXftuYOtUKpQk9dLrHHqS45JcCxwG9gC3AvdW1f1dlzuAU6ZToiSpj16BXlU/rKozgFMZPBj6SX03kGRHkoUkC4uLiyssU5I0ylifcqmqe4ErgWcCJyZZ1y06FTi4zDq7qmq+qubn5uZWVawkaXl9PuUyl+TEbvqhwAuBmxgE+8u7btuAy6dVpCRptHWju7AR2J3kOAZ/AC6pqg8nuRF4f5K/Bj4PXDjFOiVJI4wM9Kr6IvDUJdoPMDifLkk6BninqCQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEX0eQfeYJFcmuTHJDUne2LWflGRPklu61/XTL1eStJw+R+j3A39cVacDzwBel+R0YCewt6pOA/Z285KkNTIy0KvqUFVd001/m8EDok8BtgC7u267ga3TKlKSNNpY59CTbGLwfNF9wIaqOtQtuhPYMNHKJElj6R3oSR4OfBB4U1V9a3hZVRVQy6y3I8lCkoXFxcVVFStJWl6vQE/yIAZh/r6quqxrvivJxm75RuDwUutW1a6qmq+q+bm5uUnULElaQp9PuQS4ELipqv5uaNEVwLZuehtw+eTLkyT1ta5Hn2cDvw1cl+Taru3PgPOBS5JsB24HXjGdEiVJfYwM9Kr6byDLLD57suVIklbKO0UlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1Ig+37YoNW/Tzo+s2bZvO//cNdu22uIRuiQ1wkCXpEYY6JLUiD6PoHtPksNJrh9qOynJniS3dK/rp1umJGmUPkfo/wqcc1TbTmBvVZ0G7O3mJUlraGSgV9WngW8c1bwF2N1N7wa2TrguSdKYVnoOfUNVHeqm7wQ2TKgeSdIKrfqiaFUVUMstT7IjyUKShcXFxdVuTpK0jJUG+l1JNgJ0r4eX61hVu6pqvqrm5+bmVrg5SdIoKw30K4Bt3fQ24PLJlCNJWqk+H1u8CPhf4IlJ7kiyHTgfeGGSW4AXdPOSpDU08rtcquqVyyw6e8K1SJJWwTtFJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDVi5K3/x4pNOz+yJtu97fxz12S7kjQuj9AlqREGuiQ1wkCXpEYY6JLUCANdkhqxqkBPck6Sm5PsT7JzUkVJksa34kBPchzwj8CLgdOBVyY5fVKFSZLGs5oj9DOB/VV1oKq+D7wf2DKZsiRJ41pNoJ8CfG1o/o6uTZK0BqZ+p2iSHcCObva+JDev8FedDHx9MlX1lwtGdlmTunqwrvGsWV0j9jHHazzHZF25YNV1Pa5Pp9UE+kHgMUPzp3ZtP6GqdgG7VrEdAJIsVNX8an/PpFnXeKxrPNY1np/3ulZzyuVzwGlJNic5HjgPuGIyZUmSxrXiI/Squj/J64H/BI4D3lNVN0ysMknSWFZ1Dr2qPgp8dEK1jLLq0zZTYl3jsa7xWNd4fq7rSlXNYjuSpCnz1n9JasQxEeijvkIgyYOTXNwt35dk09Cyt3btNyd50YzrenOSG5N8McneJI8bWvbDJNd2PxO9WNyjrlcnWRza/u8OLduW5JbuZ9uM63rHUE1fTnLv0LKpjFeS9yQ5nOT6ZZYnyd93NX8xydOGlk1zrEbV9aqunuuSfCbJU4aW3da1X5tkYcZ1PS/JN4feqz8fWja1rwLpUdefDNV0fbc/ndQtm+Z4PSbJlV0O3JDkjUv0md0+VlVr+sPgguqtwOOB44EvAKcf1ecPgH/ups8DLu6mT+/6PxjY3P2e42ZY1/OBh3XTv3+krm7+vjUcr1cD/7DEuicBB7rX9d30+lnVdVT/NzC4kD7t8Xou8DTg+mWWvwT4GBDgGcC+aY9Vz7qedWR7DL5eY9/QstuAk9dovJ4HfHi17/+k6zqq70uBT85ovDYCT+umHwF8eYn/H2e2jx0LR+h9vkJgC7C7m74UODtJuvb3V9X3quorwP7u982krqq6sqq+281exeCz+NO2mq9ceBGwp6q+UVX3AHuAc9aorlcCF01o28uqqk8D33iALluA99bAVcCJSTYy3bEaWVdVfabbLsxu3+ozXsuZ6leBjFnXTPYtgKo6VFXXdNPfBm7ip++Yn9k+diwEep+vEPhRn6q6H/gm8Kie606zrmHbGfwVPuIhSRaSXJVk64RqGqeu3+z+eXdpkiM3gB0T49WdmtoMfHKoeVrjNcpydR9LX21x9L5VwCeSXJ3Bndiz9swkX0jysSRP7tqOifFK8jAGofjBoeaZjFcGp4KfCuw7atHM9rGfmYdEH8uS/BYwD/zqUPPjqupgkscDn0xyXVXdOqOS/gO4qKq+l+T3GPzr5qwZbbuP84BLq+qHQ21rOV7HrCTPZxDozxlqfk43Vo8G9iT5UncEOwvXMHiv7kvyEuDfgdNmtO0+Xgr8T1UNH81PfbySPJzBH5E3VdW3Jvm7x3EsHKH3+QqBH/VJsg54JHB3z3WnWRdJXgC8DXhZVX3vSHtVHexeDwCfYvCXeyZ1VdXdQ7W8G3h633WnWdeQ8zjqn8RTHK9Rlqt7mmPVS5JfYfD+bamqu4+0D43VYeBDTO4040hV9a2quq+b/ijwoCQncwyMV+eB9q2pjFeSBzEI8/dV1WVLdJndPjaNCwVjXlRYx+BiwGZ+fDHlyUf1eR0/eVH0km76yfzkRdEDTO6iaJ+6nsrgQtBpR7WvBx7cTZ8M3MKELhD1rGvj0PRvAFfVjy/CfKWrb303fdKs6ur6PYnBRarMYry637mJ5S/ynctPXrD67LTHqmddj2VwTehZR7WfADxiaPozwDkzrOsXj7x3DILxq93Y9Xr/p1VXt/yRDM6znzCr8er+298LvPMB+sxsH5vYYK9yUF7C4OrwrcDbura/YnDUC/AQ4APdDv5Z4PFD676tW+9m4MUzruu/gLuAa7ufK7r2ZwHXdTv1dcD2Gdf1N8AN3favBJ40tO7vdOO4H3jNLOvq5v8COP+o9aY2XgyO1g4BP2BwjnI78Frgtd3yMHhQy63dtudnNFaj6no3cM/QvrXQtT++G6cvdO/x22Zc1+uH9q2rGPqDs9T7P6u6uj6vZvAhieH1pj1ez2Fwjv6LQ+/VS9ZqH/NOUUlqxLFwDl2SNAEGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5Jjfh/GPgd6gOhN4MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "l = []\n",
    "for b in range(256):\n",
    "    for i in range(1000):\n",
    "        A = random.randint(0, 2**b)\n",
    "        B = random.randint(0, A)\n",
    "        l.append(len(list(lehmer_accumulate(euclid_qs(A, B)))))\n",
    "print(max(l))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(41588875349823, 1072387447758585125, 178975154770675, 4614953105182700512)]"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(lehmer_accumulate([25785, 2, 3, 1, 1, 1, 29, 23, 4, 1, 3, 3, 2, 1, 20, 1, 1, 1, 13, 2, 9, 4, 1, 1, 1, 1, 2, 3, 3, 4]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3, 2, 2]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "euclid_qs(17, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.4"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "17/5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.4"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "3 + 1 /(2 + 1/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 356,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0xc9b735c5639a1516'"
      ]
     },
     "execution_count": 356,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hex(14535145444257436950)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 358,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[12079776522, 6944202407, 5135574115, 1808628292]\n",
      "[320831736, 1018828859, 1339660595, 2358489454]\n",
      "[357461893, 1135151083, 1492612976, 2627764059]\n",
      "k = i\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 358,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lehmer_loop_3(0xe0befcaf1e3ca5f1, 0xc9b735c5639a1516, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 359,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15267531864828975732"
      ]
     },
     "execution_count": 359,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "6044159827974199924 | 2**63"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 360,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16194659139127649777"
      ]
     },
     "execution_count": 360,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "0xe0befcaf1e3ca5f1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 361,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14535145444257436950"
      ]
     },
     "execution_count": 361,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "0xc9b735c5639a1516"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 362,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[16194659139127649777, 14535145444257436950, 1659513694870212827, 1259035885295734334, 400477809574478493, 57602456572298855, 54863070140685363, 2739386431613492, 75341508415523, 27092128654664, 21157251106195, 5934877548469, 3352618460788, 2582259087681, 770359373107, 271180968360, 227997436387, 43183531973, 12079776522, 6944202407, 5135574115, 1808628292]\n",
      "[None, 1, 8, 1, 3, 6, 1, 20, 36, 2, 1, 3, 1, 1, 3, 2, 1, 5, 3, 1, 1]\n",
      "[1, 0, 1, 8, 9, 35, 219, 254, 5299, 191018, 387335, 578353, 2122394, 2700747, 4823141, 17170170, 39163481, 56333651, 320831736, 1018828859, 1339660595, 2358489454]\n",
      "[0, 1, 1, 9, 10, 39, 244, 283, 5904, 212827, 431558, 644385, 2364713, 3009098, 5373811, 19130531, 43634873, 62765404, 357461893, 1135151083, 1492612976, 2627764059]\n",
      "i = 18\n",
      "[12079776522, 6944202407, 5135574115, 1808628292]\n",
      "[320831736, 1018828859, 1339660595, 2358489454]\n",
      "[357461893, 1135151083, 1492612976, 2627764059]\n",
      "k = i + 1\n",
      "19\n",
      "1\n",
      "[16194659139127649777, 14535145444257436950, 1659513694870212827, 1259035885295734334, 400477809574478493, 57602456572298855, 54863070140685363, 2739386431613492, 75341508415523, 27092128654664, 21157251106195, 5934877548469, 3352618460788, 2582259087681, 770359373107, 271180968360, 227997436387, 43183531973, 12079776522, 6944202407, 5135574115]\n",
      "[None, 1, 8, 1, 3, 6, 1, 20, 36, 2, 1, 3, 1, 1, 3, 2, 1, 5, 3, 1]\n",
      "[1, 0, 1, 8, 9, 35, 219, 254, 5299, 191018, 387335, 578353, 2122394, 2700747, 4823141, 17170170, 39163481, 56333651, 320831736, 1018828859, 1339660595]\n",
      "[0, 1, 1, 9, 10, 39, 244, 283, 5904, 212827, 431558, 644385, 2364713, 3009098, 5373811, 19130531, 43634873, 62765404, 357461893, 1135151083, 1492612976]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(1018828859, 1339660595, 1135151083, 1492612976, False)"
      ]
     },
     "execution_count": 362,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lehmer_loop(16194659139127649777, 14535145444257436950)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 390,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[52, 1, 13, 1, 2, 4, 4, 1, 1, 1, 13, 1, 23, 1, 9, 3, 3, 5, 11, 6, 1, 3, 1, 1, 4, 2, 1, 1, 1, 1, 2, 1, 32, 1]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[(2927556694930003, 154961230633081597, 3017020641586254, 159696730135159213)]"
      ]
     },
     "execution_count": 390,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(lehmer_accumulate(euclid_qs(0x518a5cc4c55ac5b050a0831b65e827e5e39fd4515e4e094961c61509e7870814,\n",
    "                0x018a5cc4c55ac5b050a0831b65e827e5e39fd4515e4e094961c61509e7870814)[:34]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 385,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 385,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "114557320433621401 < 2927556694930003"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 391,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 391,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "5818365597666026993 > 14535145444257436950"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 392,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 392,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "10841422679839906593 > 15507080595343815048"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 393,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(37, 7, -379355176803460069, 947685836737753349)"
      ]
     },
     "execution_count": 393,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd_euclid(14535145444257436950, 5818365597666026993)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 396,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 2, 3, 11, 6, 3, 13, 1, 2, 2, 3, 4, 3, 1, 1, 3, 1, 6, 1, 4, 1, 1, 1, 1, 3, 1, 862, 632, 1, 89]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[(40154122160696118,\n",
       "  57434639988632077,\n",
       "  3613807559946635531,\n",
       "  5169026865114605016)]"
      ]
     },
     "execution_count": 396,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(lehmer_accumulate(euclid_qs(15507080595343815048, 10841422679839906593)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
